## author:    A. D?r, Robert A. Huber, Yannick Stiller
## contact:   robert.huber@sbg.ac.at
## file name: MPs_analyses_v8.R
## Context:   Understanding MP Trade Preferences
## started:   2020-04-23
## Summary:   Runs analysis for hypotheses

#. Preparation ####

#clear
rm(list=ls())

#set working dir
getwd()
# the replication_materials_LSQ.RProj should automatically set the working direction to this folder

#packages
library(tidyverse)

#load data
df <- readRDS("./data/df.rds")

df$district_year <- paste0(df$district, "-", substr(df$country_year, start = 5, stop = 8))

df_long <- df %>%
  filter(abs(merge_year - survey_year) < 3) %>%
  dplyr::select(one_of("ISO_3166_2", "id", "party_short", "gov_opp", "term_first", "country_short", "country_year", "age", "gender", "left_right", "district_clean", "district_year", "lgnic",
                       "shdi", "district_density", "msch", "state_market",
                       "educ", "inc", "hypo_USA", "hypo_EU", "hypo_Pacific", "district_magnitude"),
         contains("net"), contains("rca"), contains("agreement"), contains("share")) %>%
  gather(us_agreement:pacific_agreement, value = "ta_support", key = "agreement") %>%
  mutate(district_density = log(district_density),
         gov_opp = factor(gov_opp, levels = c("Opposition", "Government")),
         net_totl_PT_ISIC_dst = ifelse(agreement == "us_agreement", net_totl_US_ISIC_dst,
                                 ifelse(agreement == "eu_agreement", net_totl_EU_ISIC_dst,
                                        ifelse(agreement == "alba_agreement", net_totl_AL_ISIC_dst,
                                               ifelse(agreement == "pacific_agreement", net_totl_PA_ISIC_dst, NA)))),
         rca_totl_PT_ISIC_dst = ifelse(agreement == "us_agreement", rca_totl_US_ISIC_dst,
                                 ifelse(agreement == "eu_agreement", rca_totl_EU_ISIC_dst,
                                        ifelse(agreement == "alba_agreement", rca_totl_AL_ISIC_dst,
                                               ifelse(agreement == "pacific_agreement", rca_totl_PA_ISIC_dst, NA)))),
         left_right = left_right-1,
         hypo = factor(ifelse(agreement == "us_agreement", hypo_USA,
                              ifelse(agreement == "eu_agreement", hypo_EU,
                                     ifelse(agreement == "pacific_agreement", hypo_Pacific, NA))),
                       levels = c(1,0),
                       labels = c("Existing", "Hypothetical")),
         net_totl_PT_PROD_dst_endo = ifelse(agreement == "us_agreement", net_totl_US_PROD_dst_endo,
                                       ifelse(agreement == "eu_agreement", net_totl_EU_PROD_dst_endo,
                                              ifelse(agreement == "alba_agreement", net_totl_AL_PROD_dst_endo,
                                                     ifelse(agreement == "pacific_agreement", net_totl_PA_PROD_dst_endo, NA)))),
         net_totl_PT_ISIC_dst_endo = ifelse(agreement == "us_agreement", net_totl_US_ISIC_dst_endo,
                                       ifelse(agreement == "eu_agreement", net_totl_EU_ISIC_dst_endo,
                                              ifelse(agreement == "alba_agreement", net_totl_AL_ISIC_dst_endo,
                                                     ifelse(agreement == "pacific_agreement", net_totl_PA_ISIC_dst_endo, NA)))),
         rca_totl_PT_PROD_dst_endo = ifelse(agreement == "us_agreement", rca_totl_US_PROD_dst_endo,
                                       ifelse(agreement == "eu_agreement", rca_totl_EU_PROD_dst_endo,
                                              ifelse(agreement == "alba_agreement", rca_totl_AL_PROD_dst_endo,
                                                     ifelse(agreement == "pacific_agreement", rca_totl_PA_PROD_dst_endo, NA)))),
         rca_totl_PT_ISIC_dst_endo = ifelse(agreement == "us_agreement", rca_totl_US_ISIC_dst_endo,
                                       ifelse(agreement == "eu_agreement", rca_totl_EU_ISIC_dst_endo,
                                              ifelse(agreement == "alba_agreement", rca_totl_AL_ISIC_dst_endo,
                                                     ifelse(agreement == "pacific_agreement", rca_totl_PA_ISIC_dst_endo, NA)))),
         district_magnitude = ifelse(country_short == "CHL", 2, district_magnitude), # In Chile, electoral districts are nested within the regions which we use for competitiveness values. Each district elects 2 MPs.
         district_magnitude = ifelse(country_short == "ARG", district_magnitude / 2, district_magnitude), # In Argentina, only half of legislators are elected each year
         small_dist = ifelse(district_magnitude <= 5, "Small",
                             ifelse(district_magnitude > 5, "Large", NA)),
         msch_dummy = factor(ifelse(msch < 8, "Low",
                             ifelse(msch > 7, "High", NA)),
                             levels = c("Low", "High")),
         educ = ifelse(educ == "Tertiary", "Tertiary", "Non-tertiary"),
         left_right_cat = factor(car::recode(left_right, "0:1 = 'Left'; 2:3 = 'Center-left'; 
                                      4:5 = 'Center'; 6:7 = 'Center-right'; 8:9 = 'Right'"),
                                 levels = c("Left", "Center-left", "Center", "Center-right", "Right"))) %>%
  dplyr::select(-contains("_US"), -contains("_EU"), -contains("_AL"), -contains("_PA")) %>%
  filter(agreement != "alba_agreement" & !is.na(ta_support) & !grepl("-NAT", ISO_3166_2) & !grepl("-EMIGRANTS", ISO_3166_2) & !grepl("-MINORITIES", ISO_3166_2) & !grepl("-OTHER", ISO_3166_2))

# . Variables ####
dv <- "ta_support ~ "

controls <- "+ left_right + gender + agreement + country_year"
n_controls <- 4 # number of control variables + agreement dummies (without country-year dummies). Used for screenreg.

iv_total <- c("rca_totl_PT_ISIC_dst", "net_totl_PT_ISIC_dst", "rca_totl_WL_ISIC_dst", "net_totl_WL_ISIC_dst")

# . Main Models ####

df_coef <- data.frame(model = as.character(),
                      var_name = as.character(),
                      coef = as.numeric(),
                      se = as.numeric())

m_total_out <- list()
m_total <- list()
for (i in 1:length(iv_total)) {
  m_total_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls), data=df_long, model="pooling", index=c("district_year"))
  m_total[[i]] <- lmtest::coeftest(m_total_out[[i]], vcov=plm::vcovDC(m_total_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_total[i][[1]])

  df_coef_temp <- data.frame(model = iv_total[i],
                             var_name = rownames(m_total[i][[1]])[2:(2+n_controls)],
                             coef = m_total[i][[1]][2:(2+n_controls)],
                             se = m_total[i][[1]][(mod_nrow + 2):(mod_nrow + (2+n_controls))])

  df_coef <- rbind(df_coef, df_coef_temp)
}

texreg::screenreg(list(m_total_out[[1]], m_total_out[[2]], m_total_out[[3]], m_total_out[[4]]),
                  override.se = list(m_total[[1]][,2], m_total[[2]][,2], m_total[[3]][,2], m_total[[4]][,2]),
                  override.pvalues = list(m_total[[1]][,4], m_total[[2]][,4], m_total[[3]][,4], m_total[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short")
                  
texreg::texreg(list(m_total_out[[1]], m_total_out[[2]]),
               override.se = list(m_total[[1]][,2], m_total[[2]][,2]),
               override.pvalues = list(m_total[[1]][,4], m_total[[2]][,4]),
               file = "./tables/TableE1.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (Partner Competitiveness)",
               caption.above = T,
               label = "tab:main_regression",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                        "Pacific Agreement", "US Agreement",
                                        rep("Subnational Trade Competitiveness", 1)), 
               reorder.coef = c(2:6, 1))


texreg::texreg(list(m_total_out[[3]],m_total_out[[4]]),
               override.se = list(m_total[[3]][,2], m_total[[4]][,2]),
               override.pvalues = list(m_total[[3]][,4], m_total[[4]][,4]),
               file = "./tables/TableE6.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (World Competitiveness)",
               caption.above = T,
               label = "tab:app_regression_world",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the world. Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:6, 1))

df_coef %>%
  filter(grepl(x = model, pattern = "PT")) %>%
  mutate(comp_measure = if_else(grepl(x = model, pattern = "rca"), "RCA", "EX/IM"),
         agg_measure = if_else(grepl(x = model, pattern = "ISIC"), "Worker", "Product"),
         var_lab = factor(ifelse(grepl(x = var_name, pattern = "rca"), "Subnational trade\ncompetitiveness",
                                 ifelse(grepl(x = var_name, pattern = "net"), "Subnational trade\ncompetitiveness",
                                        ifelse(var_name == "left_right", "Political Ideology",
                                               ifelse(var_name == "genderFemale", "Female",
                                                      ifelse(var_name == "agreementpacific_agreement", "Pacific Agreement",
                                                             ifelse(var_name == "agreementus_agreement", "US Agreement", NA)))))),
                          levels = c("Subnational trade\ncompetitiveness", "Political Ideology", 
                                     "Female", "Pacific Agreement", "US Agreement"))) %>%
  filter(agg_measure == "Worker") %>%
  ggplot(., aes(x=forcats::fct_rev(var_lab), y = coef, colour = comp_measure, shape = comp_measure, group = paste0(agg_measure, "-", comp_measure))) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_point(position = position_dodge(width = .75), size = 2) +
  geom_linerange(aes(ymin = coef - 1.96*se, ymax = coef + 1.96*se), position = position_dodge(width = .75)) +
  geom_linerange(aes(ymin = coef - 1.645*se, ymax = coef + 1.645*se), position = position_dodge(width = .75), linewidth = 1.1) +
  ggthemes::theme_tufte() +
  coord_flip()+
  labs(x = "Variable", y = "Estimate") +
  scale_linetype_discrete("Partner:") + 
  scale_colour_discrete("Measure:") +
  scale_shape_manual("Measure:", values = c(0,19)) +
  theme(legend.position="bottom") +
  theme(text = element_text(size = 12)) +
  guides(linetype = guide_legend(nrow = 1),
         colour = guide_legend(nrow = 1),
         shape = guide_legend(nrow = 1)) +
  NULL

ggsave(filename = "./figures/Figure1.pdf",
       height = 10, 
       width = 19, units = "cm")

ggsave(filename = "./figures/Figure1.png",
       height = 10, 
       width = 19, units = "cm")

# . Robustness Checks #####
controls <- "+ left_right + gender + agreement + country_year"
controls_dst <- "+ lgnic + district_density +msch"
controls_leg <- "+ educ + inc"
controls_par <- "+ party_short"

df_coef <- data.frame(model = as.character(),
                      var_name = as.character(),
                      coef = as.numeric(),
                      se = as.numeric())

# .. Additional Controls ####

# ... District Characteristics ####

m_robust_out <- list()
m_robust <- list()
for (i in 1:length(iv_total)) {
  m_robust_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls, controls_dst), data=df_long, model="pooling", index=c("district_clean"))
  m_robust[[i]] <- lmtest::coeftest(m_robust_out[[i]], vcov=plm::vcovDC(m_robust_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_robust[i][[1]])
  
  df_coef_temp <- data.frame(model = paste0(iv_total[i], controls_dst),
                             var_name = rownames(m_robust[i][[1]])[2:3],
                             coef = m_robust[i][[1]][2:3],
                             se = m_robust[i][[1]][(mod_nrow + 2):(mod_nrow + 3)])
  
  df_coef <- rbind(df_coef, df_coef_temp)
}

texreg::screenreg(list(m_robust_out[[1]], m_robust_out[[2]], m_robust_out[[3]], m_robust_out[[4]]),
                  override.se = list(m_robust[[1]][,2], m_robust[[2]][,2], m_robust[[3]][,2], m_robust[[4]][,2]),
                  override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4], m_robust[[3]][,4], m_robust[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls+3),
                                        rep("Comp", 3)))

texreg::texreg(list(m_robust_out[[1]], m_robust_out[[2]]),
               override.se = list(m_robust[[1]][,2], m_robust[[2]][,2]),
               override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4]),
               file = "./tables/TableE7.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (District Characteristics)",
               caption.above = T,
               label = "tab:robust_dst",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement", "GNI per capita", "Log District Denisty", "District Education",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:9, 1))


# ... Legislator Characteristics ####
m_robust_out <- list()
m_robust <- list()
for (i in 1:length(iv_total)) {
  m_robust_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls, controls_leg), data=df_long, model="pooling", index=c("district_clean"))
  m_robust[[i]] <- lmtest::coeftest(m_robust_out[[i]], vcov=plm::vcovDC(m_robust_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_robust[i][[1]])
  
  df_coef_temp <- data.frame(model = paste0(iv_total[i], controls_leg),
                             var_name = rownames(m_robust[i][[1]])[2:3],
                             coef = m_robust[i][[1]][2:3],
                             se = m_robust[i][[1]][(mod_nrow + 2):(mod_nrow + 3)])
  
  df_coef <- rbind(df_coef, df_coef_temp)
}

texreg::screenreg(list(m_robust_out[[1]], m_robust_out[[2]], m_robust_out[[3]], m_robust_out[[4]]),
                  override.se = list(m_robust[[1]][,2], m_robust[[2]][,2], m_robust[[3]][,2], m_robust[[4]][,2]),
                  override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4], m_robust[[3]][,4], m_robust[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year")

texreg::texreg(list(m_robust_out[[1]], m_robust_out[[2]]),
               override.se = list(m_robust[[1]][,2], m_robust[[2]][,2]),
               override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4]),
               file = "./tables/TableE8.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (Individual Characteristics)",
               caption.above = T,
               label = "tab:robust_leg",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.map = list("(Intercept)" = NA, "rca_totl_PT_ISIC_dst" = "Subnational Trade Competitiveness", "net_totl_PT_ISIC_dst" = "Subnational Trade Competitiveness",
                                   "left_right" = "Political Ideology", "genderFemale" = "Female", "agreementpacific_agreement" = "Pacific Agreement", "agreementus_agreement" = "US Agreement",
                                   "educTertiary" = "Tertiary Education", "inc4000-7000 USD" = "Income (4-7k)", "inc7000-10000 USD" = "Income (7-10k)", "incMore than 10000 USD" = "Income (above 10k)"),
               reorder.coef = c(2:10, 1))

# ... Party Characteristics ####

df_partyeffects <- data.frame(model = as.character(),
                              coef = as.numeric())

m_robust_out <- list()
m_robust <- list()
for (i in 1:length(iv_total)) {
  m_robust_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls, controls_par), data=df_long, model="pooling", index=c("district_clean"))
  m_robust[[i]] <- lmtest::coeftest(m_robust_out[[i]], vcov=plm::vcovDC(m_robust_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_robust[i][[1]])
  
  df_coef_temp <- data.frame(model = paste0(iv_total[i], controls_par),
                             var_name = rownames(m_robust[i][[1]])[2:3],
                             coef = m_robust[i][[1]][2:3],
                             se = m_robust[i][[1]][(mod_nrow + 2):(mod_nrow + 3)])
  
  df_coef <- rbind(df_coef, df_coef_temp)
  
  df_partyeffects_temp <- data.frame(model = paste0(iv_total[i]),
                                     coef = m_robust[i][[1]][grepl(x =  rownames(m_robust[i][[1]]), "party")])
  
  df_partyeffects <- rbind(df_partyeffects, df_partyeffects_temp)
}

texreg::screenreg(list(m_robust_out[[1]], m_robust_out[[2]], m_robust_out[[3]], m_robust_out[[4]]),
                  override.se = list(m_robust[[1]][,2], m_robust[[2]][,2], m_robust[[3]][,2], m_robust[[4]][,2]),
                  override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4], m_robust[[3]][,4], m_robust[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls),
                                        rep("Comp", 3)))

texreg::texreg(list(m_robust_out[[1]], m_robust_out[[2]]),
               override.se = list(m_robust[[1]][,2], m_robust[[2]][,2]),
               override.pvalues = list(m_robust[[1]][,4], m_robust[[2]][,4]),
               file = "./tables/TableE9.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (Party Characteristics)",
               caption.above = T,
               label = "tab:robust_par",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave and party-fixed effects omitted.}",
               omit.coef = "country_year|party_short",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:6, 1))

df_coef %>%
  filter(grepl(x = model, pattern = "PT")) %>%
  mutate(comp_measure = if_else(grepl(x = model, pattern = "rca"), "RCA", "EX/IM"),
         agg_measure = if_else(grepl(x = model, pattern = "ISIC"), "Worker", "Product"),
         partner = if_else(grepl(x = model, pattern = "PT"), "Partner", "World"),
         controls = if_else(grepl(x = model, pattern = "lgnic"), "District",
                            if_else(grepl(x=model, pattern = "educ"), "Individual", "Party")),
         var_lab = factor(ifelse(grepl(x = var_name, pattern = "rca"), "Subnational trade\ncompetitiveness",
                                 ifelse(grepl(x = var_name, pattern = "net"), "Subnational trade\ncompetitiveness",
                                        ifelse(var_name == "left_right", "Political Ideology", NA))),
                          levels = c("Subnational trade\ncompetitiveness", 
                                     "Political Ideology")),
         var_name = factor(ifelse(var_lab == "Subnational trade\ncompetitiveness", paste0(comp_measure, "-", agg_measure),
                                  ifelse(var_lab == "Political Ideology", "Political Ideology", NA)),
                           levels = c("EX/IM-Product", "RCA-Product", "EX/IM-Worker", "RCA-Worker", "Political Ideology"),
                           labels = c(NA, NA, "EX/IM", "RCA", "Political Ideology"))) %>%
  filter(agg_measure == "Worker",
         var_lab != "Political Ideology") %>%
  ggplot(., aes(x=var_name, y = coef, colour = comp_measure, shape = comp_measure, group = paste0(agg_measure, "-", comp_measure))) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_point(position = position_dodge(width = .75), size = 2) +
  geom_linerange(aes(ymin = coef - 1.96*se, ymax = coef + 1.96*se), position = position_dodge(width = .75)) +
  geom_linerange(aes(ymin = coef - 1.645*se, ymax = coef + 1.645*se), position = position_dodge(width = .75), size = 1.1) +
  ggthemes::theme_tufte() +
  facet_wrap(~controls) + 
  coord_flip()+
  labs(x = "Variable", y = "Estimate") +
  scale_colour_discrete("Measure:") +
  scale_shape_manual("Measure:", values = c(0,19)) +
  theme(legend.position="none") +
  theme(text = element_text(size = 12)) +
  guides(colour = guide_legend(nrow = 1),
         shape = guide_legend(nrow = 1)) +
  NULL

ggsave(filename = "./figures/Figure4.pdf",
       height = 7,
       width = 19, units = "cm")

ggsave(filename = "./figures/Figure4.png",
       height = 7,
       width = 19, units = "cm")

# ... By agreement ####

controls_agreement <- gsub(x = paste(controls), pattern = " \\+ agreement", replacement = "")
n_controls_agreement <- n_controls - 2

# .... USA ####

df_us <- subset(df_long, agreement == "us_agreement")

df_coef_us <- data.frame(model = as.character(),
                         var_name = as.character(),
                         coef = as.numeric(),
                         se = as.numeric())

m_total_us_out <- list()
m_total_us <- list()
for (i in 1:length(iv_total)) {
  m_total_us_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls_agreement), data=df_us, model="pooling", index=c("district_clean"))
  m_total_us[[i]] <- lmtest::coeftest(m_total_us_out[[i]], vcov=plm::vcovDC(m_total_us_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_total_us[i][[1]])
  
  df_coef_temp <- data.frame(model = iv_total[i],
                             var_name = rownames(m_total_us[i][[1]])[2:(2+n_controls_agreement)],
                             coef = m_total_us[i][[1]][2:(2+n_controls_agreement)],
                             se = m_total_us[i][[1]][(mod_nrow + 2):(mod_nrow + (2+n_controls_agreement))])
  
  df_coef_us <- rbind(df_coef_us, df_coef_temp)
}

texreg::screenreg(list(m_total_us_out[[1]], m_total_us_out[[2]], m_total_us_out[[3]], m_total_us_out[[4]]),
                  override.se = list(m_total_us[[1]][,2], m_total_us[[2]][,2], m_total_us[[3]][,2], m_total_us[[4]][,2]),
                  override.pvalues = list(m_total_us[[1]][,4], m_total_us[[2]][,4], m_total_us[[3]][,4], m_total_us[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls_agreement),
                                        rep("Comp", 3)))

texreg::texreg(list(m_total_us_out[[1]], m_total_us_out[[2]]),
               override.se = list(m_total_us[[1]][,2], m_total_us[[2]][,2]),
               override.pvalues = list(m_total_us[[1]][,4], m_total_us[[2]][,4]),
               file = "./tables/TableE10.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (US Agreement)",
               caption.above = T,
               label = "tab:agreement_us",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year|party_short",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:4,1))

# .... EU ####

df_eu <- subset(df_long, agreement == "eu_agreement")

df_coef_eu <- data.frame(model = as.character(),
                         var_name = as.character(),
                         coef = as.numeric(),
                         se = as.numeric())

m_total_eu_out <- list()
m_total_eu <- list()
for (i in 1:length(iv_total)) {
  m_total_eu_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls_agreement), data=df_eu, model="pooling", index=c("district_year"))
  m_total_eu[[i]] <- lmtest::coeftest(m_total_eu_out[[i]], vcov=plm::vcovDC(m_total_eu_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_total_eu[i][[1]])
  
  df_coef_temp <- data.frame(model = iv_total[i],
                             var_name = rownames(m_total_eu[i][[1]])[2:(2+n_controls_agreement)],
                             coef = m_total_eu[i][[1]][2:(2+n_controls_agreement)],
                             se = m_total_eu[i][[1]][(mod_nrow + 2):(mod_nrow + (2+n_controls_agreement))])
  
  df_coef_eu <- rbind(df_coef_eu, df_coef_temp)
}

texreg::screenreg(list(m_total_eu_out[[1]], m_total_eu_out[[2]], m_total_eu_out[[3]], m_total_eu_out[[4]]),
                  override.se = list(m_total_eu[[1]][,2], m_total_eu[[2]][,2], m_total_eu[[3]][,2], m_total_eu[[4]][,2]),
                  override.pvalues = list(m_total_eu[[1]][,4], m_total_eu[[2]][,4], m_total_eu[[3]][,4], m_total_eu[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls_agreement),
                                        rep("Comp", 3)))

texreg::texreg(list(m_total_eu_out[[1]], m_total_eu_out[[2]]),
               override.se = list(m_total_eu[[1]][,2], m_total_eu[[2]][,2]),
               override.pvalues = list(m_total_eu[[1]][,4], m_total_eu[[2]][,4]),
               file = "./tables/TableE11.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (EU Agreement)",
               caption.above = T,
               label = "tab:agreement_eu",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year|party_short",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:4,1))

# .... Pacific ####

df_pa <- subset(df_long, agreement == "pacific_agreement")

df_coef_pa <- data.frame(model = as.character(),
                         var_name = as.character(),
                         coef = as.numeric(),
                         se = as.numeric())

m_total_pa_out <- list()
m_total_pa <- list()
for (i in 1:length(iv_total)) {
  m_total_pa_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls_agreement), data=df_pa, model="pooling", index=c("district_year"))
  m_total_pa[[i]] <- lmtest::coeftest(m_total_pa_out[[i]], vcov=plm::vcovDC(m_total_pa_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_total_pa[i][[1]])
  
  df_coef_temp <- data.frame(model = iv_total[i],
                             var_name = rownames(m_total_pa[i][[1]])[2:(2+n_controls_agreement)],
                             coef = m_total_pa[i][[1]][2:(2+n_controls_agreement)],
                             se = m_total_pa[i][[1]][(mod_nrow + 2):(mod_nrow + (2+n_controls_agreement))])
  
  df_coef_pa <- rbind(df_coef_pa, df_coef_temp)
}

texreg::screenreg(list(m_total_pa_out[[1]], m_total_pa_out[[2]], m_total_pa_out[[3]], m_total_pa_out[[4]]),
                  override.se = list(m_total_pa[[1]][,2], m_total_pa[[2]][,2], m_total_pa[[3]][,2], m_total_pa[[4]][,2]),
                  override.pvalues = list(m_total_pa[[1]][,4], m_total_pa[[2]][,4], m_total_pa[[3]][,4], m_total_pa[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls_agreement),
                                        rep("Comp", 3)))

texreg::texreg(list(m_total_pa_out[[1]], m_total_pa_out[[2]]),
               override.se = list(m_total_pa[[1]][,2], m_total_pa[[2]][,2]),
               override.pvalues = list(m_total_pa[[1]][,4], m_total_pa[[2]][,4]),
               file = "./tables/TableE12.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (Pacific Alliance)",
               caption.above = T,
               label = "tab:agreement_pa",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year|party_short",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     rep("Subnational Trade Competitiveness", 1)),
               reorder.coef = c(2:4,1))

# .. Non-linearity of interaction terms ####

df_long$us_agreement <- factor(ifelse(df_long$agreement == "us_agreement", 1, 0))
df_long$eu_agreement <- factor(ifelse(df_long$agreement == "eu_agreement", 1, 0))

library(interflex)

intdata <- df_long[,c("ta_support",
                      "rca_totl_PT_ISIC_dst",
                      "left_right",
                      "country_year",
                      "gender",
                      "us_agreement",
                      "eu_agreement", 
                      "district_magnitude",
                      "district_year")] %>% 
  mutate(ldm = log(district_magnitude)) %>% 
  select(-district_magnitude)

intdata <- na.omit(intdata)
intdata <- as.data.frame(intdata)

# ... RCA ####

s5.kernela <-interflex(Y = "ta_support", 
                       D = "rca_totl_PT_ISIC_dst", 
                       X = "left_right",
                       FE = "country_year",
                       Z = c("gender","us_agreement","eu_agreement", "ldm"),
                       data = intdata, 
                       estimator = "binning",
                       vcov.type = "cluster",
                       cl = "district_year",
                       nbins = 5,
                       base = "Stable",
                       Ylabel = "Trade support\n",Dlabel = "RCA",
                       Xlabel = "Left-Right Self-Placement",
                       show.grid = F,theme.bw = T)

s5.kernela$figure

ggsave(filename = "./figures/FigureE2.pdf",
       height = 12, 
       width = 19, units = "cm")

s5.kernelb <-interflex(Y = "ta_support", 
                       D = "rca_totl_PT_ISIC_dst", 
                       X = "left_right",
                       FE = "country_year",
                       Z = c("gender","us_agreement","eu_agreement", "ldm"),
                       data = intdata, 
                       estimator = "kernel",
                       bw = 1,
                       vcov.type = "cluster",
                       cl = "district_year",
                       base = "Stable",Ylabel = "Trade support\n",Dlabel = "RCA",
                       Xlabel = "Left-Right Self-Placement",
                       show.grid = F,theme.bw = T)

s5.kernelb$figure

ggsave(filename = "./figures/FigureE3.pdf",
       height = 12, 
       width = 19, units = "cm")

s5.kernelc <-interflex(Y = "ta_support", 
                       D = "rca_totl_PT_ISIC_dst", 
                       X = "ldm",
                       FE = "country_year",
                       Z = c("gender","us_agreement","eu_agreement", "left_right"),
                       data = intdata, 
                       estimator = "binning",
                       vcov.type = "cluster",
                       cl = "district_year",
                       nbins = 5,
                       base = "Stable",Ylabel = "Trade support\n",Dlabel = "RCA",
                       Xlabel = "Left-Right Self-Placement",
                       show.grid = F,theme.bw = T)

s5.kernelc$figure

ggsave(filename = "./figures/FigureE4.pdf",
       height = 12, 
       width = 19, units = "cm")

s5.kerneld <-interflex(Y = "ta_support", 
                       D = "rca_totl_PT_ISIC_dst",  
                       X = "ldm",
                       FE = "country_year",
                       Z = c("gender","us_agreement","eu_agreement", "left_right"),
                       data = intdata, 
                       estimator = "kernel",
                       bw = 1,
                       vcov.type = "cluster",
                       cl = "district_year",
                       base = "Stable",Ylabel = "Trade support\n",Dlabel = "RCA",
                       Xlabel = "Left-Right Self-Placement",
                       show.grid = F,theme.bw = T)

s5.kerneld$figure

ggsave(filename = "./figures/FigureE5.pdf",
       height = 12, 
       width = 19, units = "cm")

# . Exploration ####

# .. Interactions ####

# ... District magnitude & ideology ####

df_marg_size <- data.frame(iv_name = as.character(),
                           factor = as.character(),
                           district_magnitude = as.numeric(),
                           AME = as.numeric(),
                           SE = as.numeric(),
                           p = as.numeric())

df_marg_ideo <- data.frame(iv_name = as.character(),
                           factor = as.character(),
                           left_right_cat = as.numeric(),
                           AME = as.numeric(),
                           SE = as.numeric(),
                           p = as.numeric())

m_comb_total_lm <- list()
m_comb_total_test <- list()
for (i in 1:length(iv_total)) {
  df_long$comp_temp <- df_long[,iv_total[i]]
  
  m_comb_total_lm[[i]] <- lm(formula=paste0(dv, "comp_temp * log(district_magnitude) + comp_temp*I(log(district_magnitude)^2) + comp_temp * left_right + gender + agreement + country_year"), data=df_long)
  m_comb_total <- plm::plm(formula=paste0(dv, "comp_temp * log(district_magnitude) + comp_temp*I(log(district_magnitude)^2) + comp_temp * left_right + gender + agreement + country_year"), data=df_long, model="pooling", index=c("district_clean"))
  m_comb_total_test[[i]] <- lmtest::coeftest(m_comb_total, vcov=plm::vcovDC(m_comb_total, type="HC0"))
  
  margins_size <- margins::margins(model = m_comb_total_lm[[i]], data = df_long,
                                   vcov = plm::vcovDC(m_comb_total, type="HC0"),
                                   at = list(district_magnitude = c(1:20)),
                                   variables = c("comp_temp"))
  
  margins_ideo <- margins::margins(model = m_comb_total_lm[[i]], data = df_long,
                                   vcov = plm::vcovDC(m_comb_total, type="HC0"),
                                   at = list(left_right = c(0:9)),
                                   variables = c("comp_temp"))
  
  df_marg_temp_size <- data.frame(iv_name = iv_total[i],
                                  factor = rep(summary(margins_size)["factor"]),
                                  district_magnitude = c(summary(margins_size)["district_magnitude"]),
                                  AME = c(summary(margins_size)["AME"]),
                                  SE = c(summary(margins_size)["SE"]),
                                  p = c(round(summary(margins_size)["p"], 3)))
  
  df_marg_size <- rbind(df_marg_size, df_marg_temp_size)
  
  df_marg_temp_ideo <- data.frame(iv_name = iv_total[i],
                                  factor = rep(summary(margins_ideo)["factor"]),
                                  left_right = c(summary(margins_ideo)["left_right"]),
                                  AME = c(summary(margins_ideo)["AME"]),
                                  SE = c(summary(margins_ideo)["SE"]),
                                  p = c(round(summary(margins_ideo)["p"], 3)))
  
  df_marg_ideo <- rbind(df_marg_ideo, df_marg_temp_ideo)
  
  df_long$comp_temp <- NULL
}

df_marg_size %>%
  mutate(comp_measure = factor(if_else(grepl(x = iv_name, pattern = "rca"), "RCA", "EX/IM"),
                               levels = c("RCA", "EX/IM")),
         agg_measure = if_else(grepl(x = iv_name, pattern = "ISIC"), "Worker", "Product"),
         partner = if_else(grepl(x = iv_name, pattern = "PT"), "Partner", "World"),
         group = paste0(agg_measure,"\n",partner),
         factor = factor(factor,
                         levels = c("comp_temp"),
                         labels = c("Competitiveness"))) %>%
  filter(partner == "Partner" & agg_measure == "Worker") %>%
  ggplot(., aes(x = district_magnitude, y = AME)) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_point(size = 2) +
  geom_linerange(aes(ymin = AME - 1.96*SE, ymax = AME + 1.96*SE)) +
  geom_linerange(aes(ymin = AME - 1.645*SE, ymax = AME + 1.645*SE), size = 1.1) +
  facet_grid(~comp_measure) + 
  ggthemes::theme_tufte() + 
  labs(x="District Magnitude", y="Marginal Effect") +
  theme(legend.position = "bottom") +
  theme(text = element_text(size = 12)) +
  scale_x_continuous(breaks = c(1:20)) +
  NULL

ggsave(filename = "./figures/Figure2.pdf",
       height = 7, 
       width = 19, units = "cm")

ggsave(filename = "./figures/Figure2.png",
       height = 7, 
       width = 19, units = "cm")

df_marg_ideo %>%
  mutate(comp_measure = factor(if_else(grepl(x = iv_name, pattern = "rca"), "RCA", "EX/IM"),
                               levels = c("RCA", "EX/IM")),
         agg_measure = if_else(grepl(x = iv_name, pattern = "ISIC"), "Worker", "Product"),
         partner = if_else(grepl(x = iv_name, pattern = "PT"), "Partner", "World"),
         group = paste0(agg_measure,"\n",partner),
         factor = factor(factor,
                         levels = c("comp_temp", "left_right"),
                         labels = c("Competitiveness", "Political Ideology"))) %>%
  filter(partner == "Partner" & agg_measure == "Worker",
         factor != "Political Ideology") %>%
  ggplot(., aes(x = left_right, y = AME)) + 
  geom_hline(yintercept = 0, lty = "dashed") + 
  geom_point(position = position_dodge(width = .75), size = 2) +
  geom_linerange(aes(ymin = AME - 1.96*SE, ymax = AME + 1.96*SE), position = position_dodge(width = .75)) +
  geom_linerange(aes(ymin = AME - 1.645*SE, ymax = AME + 1.645*SE), position = position_dodge(width = .75), size = 1.1) +
  facet_grid(~comp_measure) + 
  ggthemes::theme_tufte() + 
  labs(x="Political Ideology", y="Marginal Effect") +
  scale_x_continuous(breaks = 0:9) +
  scale_shape_manual("Political Ideology:", values = c(0,19)) +
  scale_colour_discrete("Political Ideology:") + 
  theme(legend.position = "none") +
  theme(text = element_text(size = 12)) +
  NULL

ggsave(filename = "./figures/Figure3.pdf",
       height = 7, 
       width = 19, units = "cm")

ggsave(filename = "./figures/Figure3.png",
       height = 7, 
       width = 19, units = "cm")

texreg::texreg(list(m_comb_total_lm[[1]], m_comb_total_lm[[2]]),
               override.se = list(m_comb_total_test[[1]][,2], m_comb_total_test[[2]][,2]),
               override.pvalues = list(m_comb_total_test[[1]][,4], m_comb_total_test[[2]][,4]),
               file = "./tables/TableE2.tex",
               caption = "Subnational Trade Competitiveness, Boundary Conditions, and Trade Attitudes",
               caption.above = T,
               label = "tab:exploration_comb",
               stars = c(.1, .05, .01),
               scalebox = 0.9,
               dcolumn = T,
               use.package = F,
               #sideways = T,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               reorder.coef = c(2:5, 9:11, 6:8, 1),
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Log. District Magnitude", "Log. District Magnitude Sq.",
                                     "Political Ideology",  "Female",
                                     "Pacific Agreement", "US Agreement", "Comp. x Log. Dst. Magnitude", "Comp. x Dst. Log. Magnitude Sq.", "Comp. x Pol. Ideology"))

# .. Open/Closed Lists ####

df_long <- dplyr::left_join(df_long, data.frame(country_short = unique(df$country_short),
                                                list_type = c("Closed", "Closed", "Open",
                                                              "Closed", "Closed", "Closed",
                                                              "Open", "Closed", "Open",
                                                              "Closed", "Closed", "Closed",
                                                              "Open", "Closed", "Closed",
                                                              "Closed", "Closed")))

df_marg <- data.frame(iv_name = as.character(),
                      lgnic = as.numeric(),
                      AME = as.numeric(),
                      SE = as.numeric(),
                      p = as.numeric())

m_list_total_lm <- list()
m_list_total_test <- list()
#plot_total <- list()
for (i in 1:length(iv_total)) {
  df_long$comp_temp <- df_long[,iv_total[i]]
  
  m_list_total_lm[[i]] <- lm(formula=paste0(dv, "comp_temp", " * list_type", controls), data=df_long)
  m_list_total <- plm::plm(formula=paste0(dv, "comp_temp", " * list_type", controls), data=df_long, model="pooling", index=c("district_year"))
  m_list_total_test[[i]] <- lmtest::coeftest(m_list_total, vcov=plm::vcovDC(m_list_total, type="HC0"))
  
  margins <- margins::margins(model = m_list_total_lm[[i]], data = df_long,
                              vcov = plm::vcovDC(m_list_total, type="HC0"),
                              at = list(list_type = c("Closed", "Open")),
                              variables = "comp_temp")
  
  df_marg_temp <- data.frame(iv_name = rep(iv_total[i],2),
                             hypo = c(summary(margins)["list_type"]),
                             AME = c(summary(margins)["AME"]),
                             SE = c(summary(margins)["SE"]),
                             p = c(round(summary(margins)["p"], 3)))
  
  df_marg <- rbind(df_marg, df_marg_temp)
  
  df_long$comp_temp <- NULL
}

texreg::screenreg(list(m_list_total_lm[[1]], m_list_total_lm[[2]], m_list_total_lm[[3]], m_list_total_lm[[4]]),
                  override.se = list(m_list_total_test[[1]][,2], m_list_total_test[[2]][,2], m_list_total_test[[3]][,2], m_list_total_test[[4]][,2]),
                  override.pvalues = list(m_list_total_test[[1]][,4], m_list_total_test[[2]][,4], m_list_total_test[[3]][,4], m_list_total_test[[4]][,4]),
                  stars = c(.1, .05, .01),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls+1), "Comp x Hypo"))

texreg::texreg(list(m_list_total_lm[[1]],m_list_total_lm[[2]]),
               override.se = list(m_list_total_test[[1]][,2], m_list_total_test[[2]][,2]),
               override.pvalues = list(m_list_total_test[[1]][,4], m_list_total_test[[2]][,4]),
               file = "./tables/TableE13.tex",
               caption = "Subnational Trade Competitiveness, Open List, and Trade Attitudes",
               caption.above = T,
               label = "tab:list_type",
               stars = c(.1, .05, .01),
               scalebox = 0.9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Open List", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement", "Comp. x Open List"),
               reorder.coef = c(2:3,8, 4:7, 1))

# . Endogeneity ####

# .. Hypothetical Argeements ####

df_pred <- data.frame(iv_name = as.character(),
                      comp = as.numeric(),
                      gov_opp = as.numeric(),
                      pred = as.numeric(),
                      se = as.numeric())

df_marg <- data.frame(iv_name = as.character(),
                      lgnic = as.numeric(),
                      AME = as.numeric(),
                      SE = as.numeric(),
                      p = as.numeric())

m_hypo_total_lm <- list()
m_hypo_total_test <- list()

for (i in 1:length(iv_total)) {
  df_long$comp_temp <- df_long[,iv_total[i]]
  
  m_hypo_total_lm[[i]] <- lm(formula=paste0(dv, "comp_temp", " * hypo", controls), data=df_long)
  m_hypo_total <- plm::plm(formula=paste0(dv, "comp_temp", " * hypo", controls), data=df_long, model="pooling", index=c("district_year"))
  m_hypo_total_test[[i]] <- lmtest::coeftest(m_hypo_total, vcov=plm::vcovDC(m_hypo_total, type="HC0"))
  
  pred <- prediction::prediction_summary(model = m_hypo_total_lm[[i]],
                                         vcov = plm::vcovDC(m_hypo_total, type="HC0"),
                                         at = list(comp_temp = seq(from = 0,
                                                                   to = 1,
                                                                   length.out = 11),
                                                   hypo = c("Existing", "Hypothetical")),
                                         calculate_se = TRUE)
  
  df_pred_temp <- data.frame(iv_name = iv_total[i],
                             comp = pred$`at(comp_temp)`,
                             hypo = pred$`at(hypo)`,
                             pred = pred$Prediction,
                             se = pred$SE)
  
  df_pred <- rbind(df_pred, df_pred_temp)
  
  margins <- margins::margins(model = m_hypo_total_lm[[i]], data = df_long,
                              vcov = plm::vcovDC(m_hypo_total, type="HC0"),
                              at = list(hypo = c("Existing", "Hypothetical")),
                              variables = "comp_temp")
  
  df_marg_temp <- data.frame(iv_name = rep(iv_total[i],2),
                             hypo = c(summary(margins)["hypo"]),
                             AME = c(summary(margins)["AME"]),
                             SE = c(summary(margins)["SE"]),
                             p = c(round(summary(margins)["p"], 3)))
  
  df_marg <- rbind(df_marg, df_marg_temp)
  
  df_long$comp_temp <- NULL
}

texreg::screenreg(list(m_hypo_total_lm[[1]], m_hypo_total_lm[[2]], m_hypo_total_lm[[3]], m_hypo_total_lm[[4]]),
                  override.se = list(m_hypo_total_test[[1]][,2], m_hypo_total_test[[2]][,2], m_hypo_total_test[[3]][,2], m_hypo_total_test[[4]][,2]),
                  override.pvalues = list(m_hypo_total_test[[1]][,4], m_hypo_total_test[[2]][,4], m_hypo_total_test[[3]][,4], m_hypo_total_test[[4]][,4]),
                  stars = c(.1, .05, .01),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"), 
                  omit.coef = "country_year",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls+1), "Comp x Hypo"))

texreg::texreg(list(m_hypo_total_lm[[1]],m_hypo_total_lm[[2]]),
               override.se = list(m_hypo_total_test[[1]][,2], m_hypo_total_test[[2]][,2]),
               override.pvalues = list(m_hypo_total_test[[1]][,4], m_hypo_total_test[[2]][,4]),
               file = "./tables/TableE3.tex",
               caption = "Subnational Trade Competitiveness, Trade Attitudes and Agreement Status",
               caption.above = T,
               label = "tab:endo_hypo",
               stars = c(.1, .05, .01),
               scalebox = 0.9,
               dcolumn = T,
               use.package = F,
               #sideways = T,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Hypothetical", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement", "Comp. x Hypothetical"),
               reorder.coef = c(2,3,8, 4:7, 1))

# .. Only First Times ####

m_first_total_out <- list()
m_first_total <- list()
for (i in 1:length(iv_total)) {
  m_first_total_out[[i]] <- plm::plm(formula=paste0(dv, iv_total[i], controls), data=df_long, subset = term_first == "Yes", model="pooling", index=c("district_year"))
  m_first_total[[i]] <- lmtest::coeftest(m_first_total_out[[i]], vcov=plm::vcovDC(m_first_total_out[[i]], type="HC0"))
}

texreg::screenreg(list(m_first_total_out[[1]], m_first_total_out[[2]], m_first_total_out[[3]], m_first_total_out[[4]]),
                  override.se = list(m_first_total[[1]][,2], m_first_total[[2]][,2], m_first_total[[3]][,2], m_first_total[[4]][,2]),
                  override.pvalues = list(m_first_total[[1]][,4], m_first_total[[2]][,4], m_first_total[[3]][,4], m_first_total[[4]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "EX/IM_ISIC_PRT", "RCA_ISIC_WLD", "EX/IM_ISIC_WLD"),
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "Comp", rep(NA, n_controls),
                                        rep("Comp", 3)))

texreg::texreg(list(m_first_total_out[[1]], m_first_total_out[[2]]),
               override.se = list(m_first_total[[1]][,2], m_first_total[[2]][,2]),
               override.pvalues = list(m_first_total[[1]][,4], m_first_total[[2]][,4]),
               file = "./tables/TableE4.tex",
               caption = "Subnational Trade Competitiveness and Trade Attitudes (Only First Term)",
               caption.above = T,
               label = "tab:endo_first",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement", "Subnational Trade Competitiveness"),
               reorder.coef = c(2,3:6,1))


# .. Trade Attitudes on Competitiveness ####

dv_total <- paste0(iv_total, "_endo")[1:2] #we only run this analysis for the partner operationalisation

df_coef <- data.frame(model = as.character(),
                      var_name = as.character(),
                      coef = as.numeric(),
                      se = as.numeric())

m_total_out <- list()
m_total <- list()
for (i in 1:2) {
  m_total_out[[i]] <- plm::plm(formula=paste0(dv_total[i],"~ ta_support + ", iv_total[i], controls), data=df_long, model="pooling", index=c("district_year"))
  m_total[[i]] <- lmtest::coeftest(m_total_out[[i]], vcov=plm::vcovDC(m_total_out[[i]], type="HC0"))
  
  mod_nrow <- nrow(m_total[i][[1]])
  
  df_coef_temp <- data.frame(model = iv_total[i],
                             var_name = rownames(m_total[i][[1]])[2:(2+n_controls)],
                             coef = m_total[i][[1]][2:(2+n_controls)],
                             se = m_total[i][[1]][(mod_nrow + 2):(mod_nrow + (2+n_controls))])
  
  df_coef <- rbind(df_coef, df_coef_temp)
}

texreg::screenreg(list(m_total_out[[1]], m_total_out[[2]]),
                  override.se = list(m_total[[1]][,2], m_total[[2]][,2]),
                  override.pvalues = list(m_total[[1]][,4], m_total[[2]][,4]),
                  stars = c(.1, .05),
                  custom.model.names = c("RCA_ISIC_PRT", "RCA_PROD_PRT"), 
                  omit.coef = "country_year|party_short",
                  custom.coef.names = c(NA, "TA_Support", "Comp", rep(NA, n_controls),
                                        rep("Comp", 1)))

texreg::texreg(list(m_total_out[[1]], m_total_out[[2]]),
               override.se = list(m_total[[1]][,2], m_total[[2]][,2]),
               override.pvalues = list(m_total[[1]][,4], m_total[[2]][,4]),
               file = "./tables/TableE5.tex",
               caption = "Endogeneity Tests",
               caption.above = T,
               label = "tab:endo_t2",
               stars = c(.1, .05, 0.01),
               scalebox = .9,
               dcolumn = T,
               use.package = F,
               float.pos = "htb",
               #sideways = T,
               center = T,
               custom.model.names = c("RCA", "EX/IM"),
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on the district-year. RCA or EX/IM refers to the revealed comparative advantage measure and trade balance (net-trade) measure of subnational trade competitiveness. RCA or EX/IM measures are calculated vis-a-vis the respective partner (US, EU or Pacific Alliance). Country-wave fixed effects omitted.}",
               omit.coef = "country_year",
               custom.coef.names = c(NA, "Trade Support", "Subnational Trade Competitiveness", "Political Ideology", "Female",
                                     "Pacific Agreement", "US Agreement",
                                     rep("Subnational Trade Competitiveness", 1)))

